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Abstract 

Parametric imaging of the cerebral metabolic rate for glucose (CMRGlc) 
using [^^F]-fluorodeoxyglucose positron emission tomography is considered. 
Traditional imaging is hindered due to low signal to noise ratios at individual 
voxels. We propose to minimize the total variation of the tracer uptake rates 
while requiring good fit of traditional Patlak equations. This minimization 
guarantees spatial homogeneity within brain regions and good distinction 
between brain regions. Brain phantom simulations demonstrate significant 
improvement in quality of images by the proposed method as compared to 
Patlak images with post-filtering using Gaussian or median filters. 

Key words: Total variation; graphical analysis; Patlak plot; PET 
quantification; Parametric imaging; FDG; Alzheimer's disease; uptake rate. 



1. Introduction 

We focus on positron emission tomography (PET) parametric imaging 
for estimating the cerebral metabolic rate of glucose (CMRGlc) using the 
[^'^Fj-fluorodeoxyglucose (FDG) tracer. The ability to derive accurate pa- 
rameters depends upon the quality of data, the quantification method and 
the numerical algorithm. In this study, we refer to the time activity curve 
(TAG) from a given tissue location as the output, tissue TAG or TTAG, 
and the TAG from the blood pool (image-derived or arterial blood-sampled) 
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as the input, plasma TAG or PTAC. Most existing quantification meth- 
ods perform well for regions of interest (ROIs), but are not good for voxel 
level quantification due to the high level of noise. These include graphical 
methods, 20, 18], linear least squares, the weighted integration method, [1], 
generalized linear least squares, 0, 0], nonlinear least squares (NLS) and 
weighted NLS. 

All the algorithms listed above perform the quantification at each voxel 
location separately; they do not consider the kinetic similarities among 
neighboring voxels within functionally-defined regions. Thus, voxel-by- voxel 
variation in a functionally-homogeneous region may be large because of noise 
in the data. But, by incorporating the spatial constraint that parameters 
in a functionally-homogeneous region should be similar, in any of the above 
methods, the quality of the resulting parametric image for the CMRGlc 



may be improved. Zhou et al, [25[, for example, improved the parametric 



image quality by ridge regression with constraints on the rate constants. 
There, the estimation of parameters uses a linear components decomposi- 
tion of the kinetics in which each component represents a functional kinetic 
curve generated by clustering the TTACs, and the problem is solved voxel 
by voxel. Here we propose the use of a total variation (TV) penalty term 
which imposes spatial consistency between neighboring voxels. 

The total variation penalty was first introduced in the context of image 
deblurring by Rudin et al, [22]. TV can significantly suppress noise while 
recovering sharp edges because it does not penalize discontinuities. It has 
received much theoretical research attention and been utilized in many signal 
and image processing applications. While it was introduced for PET image 
reconstruction by Jonsson et al, 13], and Kisilev et al, [iH], it has apparently 
not been applied for parametric PET imaging. Instead of calculating the 
uptake rate for each voxel by Patlak's method, we propose to minimize the 
TV of the uptake rate over the entire image while also maintaining a good 
least squares fit for the Patlak equations at all voxels. Thus the parameters 
of the whole image are spatially related by the TV and solved simultaneously. 
The resultant parametric image is expected to have spatial homogeneity over 
brain regions with similar kinetics and distinct edges between brain regions 
that have different kinetics. This is validated by phantom simulations. 

In addition to proposing the new model with the TV penalty, we also pay 
careful attention to the computational complexity of the algorithm by taking 
advantage of implementations for large scale sparse matrix computations. In 
contrast to approximating the Hessian matrix, as is typical for quasi-Newton 
methods, our algorithm explicitly and accurately calculates both gradient 
and Hessian terms. The Hessian is efficiently recalculated at each iteration 
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because of its sparsity. The Quasi-Minimal Residual (QMR) method, [71], 
is used to solve the resulting large-scale linear systems. With this efficient 
implementation, the procedure described here is computationally feasible. 

The rest of the paper is organized as follows: The new algorithm with 
TV penalty is introduced in section [21 with relevant computational issues 
detailed in the appendix. The experimental data sets are described in sec- 
tion [3] and results reported in section [H Issues relevant to the proposed 
TV-Patlak method and computational aspects are discussed in section [5l 
Conclusions are presented in section [6l 



2. Methods 

The Patlak plot has been developed for systems with irreversible trap- 
ping [i^]. Most often it is applied for the analysis of FDG. The measured 
TTAC undergoes a transformation and is plotted against a normalized time. 
It is given by the expression 

gr(t)^^ /o-Cp(.)d. 
Cp{t) Cp{t) 

where Cicit) is the measured TTAC, (in counts/min/g) and Cp is the PTAC 
(in counts/min/ml) , i.e. the FDG concentration in plasma. For systems 
with irreversible compartments this plot yields a straight line after suffi- 
cient equilibration time. For the FDG tracer, the slope K represents the 
uptake rate which, together with lumped constant (LC) and glucose concen- 
tration in plasma (Cpg) allows easy calculation of the CMRGlc=ErCpg/LC, 
(in mg/min/lOOg). The intercept V is given by Vq + vB where Vq is the 
distribution volume of the reversible compartment and vB is the fractional 
blood volume. 

The linear relationship ([1]) can be rewritten as 



J^Cp{s)ds^K + Cp{t)V « CT{t), 



(2) 



and, assuming m dynamic frames over the period of equilibration, its dis- 
cretized version is 







Cp{s)dsjK + Cp{tj)V ^CT{tj), j = l,---,m. (3) 
In matrix format, 

^(^)-b, (4) 
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where ^ is a m-by-2 matrix, 



£Cp{s)ds, Cp(i2) 



V /J'"Cp(s)ds, Cp(t^) / 

and vector b = (CT(ti), CT(t2); • ' ' ■,CT{'tm))'^ ■ If we were to solve (j3|) for 
each voxel independently we would obtain a parametric image lacking spa- 
tial homogeneity and with low signal-to-noise ratio (SNR). Image denoising 
techniques could then be applied as a separate task to improve the image 
quality. Instead, obtaining all voxel parameters as a result of a global op- 
timization algorithm with a TV penalty for the entire image, the necessity 
for postprocessing should be eliminated. 

Limiting the discussion here to 2D images (although our application of 
TV is 3D), we select the active voxels to be quantified by the application of 
a brain mask, yielding a total of N voxels. Equation ^ holds with common 
matrix A dependent on Cp{t) for each voxel i, but with K, V and b replaced 
by K^''\ V^^\ and bj, respectively, where bj is obtained from the TTAC for 
voxel i. Collecting the unknowns of these voxels in vectors of uptake rates 
and intercepts 

x = (i^«,i^(2),... ,KWf and y = (y(i),y(2),... 

and requiring (j4|) in the least squares sense over all voxels, while maintaining 
minimal TV of the uptake rate for the selected image voxels, yields the global 
minimization problem 

(TV-Patlak): min$(x;y), (5) 

N 

$(x;y) = M^y + a \\Wi{A{x,, -hi)\\l (6) 

i=l 

Here the total variation norm is given by ||x||tv,/3 = '^iLi </'i(x) with (^j(x) = 
1/ (xi — Xi^Y + {xi — Xi^y + and Xi^ and Xi^ are the values associated 
with voxels to the right and below voxel i. Theoretically the TV norm is 



||x||tv,05 which is a seminorm on a space of bounded variation, 2J]. The 
small constant P is used to avoid the numerical difficulty due to the lack of 
differentiability at the origin of (pi for /? = 0. The diagonal weight matrix 
for voxel i is given by 



Wi = diag . / ^ ], j = !,■■■ ,m, (7) 
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where Aij is the scan duration of the frame at time tj, A is the tracer's decay 
constant and bij is the value of the i^^ TTAC at frame j. This weighting is 
consistent with using a simulation with variance 

Var(CT(t,)) = 5c^I^^, (8) 



13]. Sc is a common scale factor that need not be made explicit here because 
it is absorbed into the parameter a in ([6]). 

The objective function in (EJis convex, and can be solved using a stan- 
dard Newton-type algorithm, [2j] Chapter 8. To simplify the expressions 
we introduce the vector z = [x;y]. 

Algorithm 1. Given initial guess z = [x;y] and tolerance e > 
Repeat 

1. Solve for Az in 

V2$(z)Az = -V«>(z). (9) 

2. Stop if |V^>^Az| < e. 

3. Line search: Choose step size s. 

4. Update: z = z + sAz. 

Further details on the calculation of gradient vector V$(z) and Hessian 
matrix V^$(z) are provided in the appendix. Some other aspects of the 
algorithm, including discussion about constants a and /3, here chosen to be 
0.2 and 10~^, respectively, as well as other approaches for the solution of 
the TV problem are discussed in section [5l 

3. Experimental Data 

To validate the proposed parametric imaging method we performed ex- 
periments with simulated data. The MRI-based high-resolution Zubal head 



phantom is used to define the brain structures, [2a]. Each voxel in the slice 
of 256 X 256 voxels is of size 1.5 x 1.5mm^. There are 128 slices for the entire 
head and 62 defined anatomical, neurological, and taxonomical structures. 
The 11 regions for slice 64, representing a total 13708 voxels, are given in 
table [TJ For the purposes of the simulations, well-accepted values of the ki- 
netic rate parameters of these structures for the two-tissue compartmental 



model of the FDG tracer [Ij], are assigned. Notice that in order to better 
model the true biochemical process we assume k/^ > 0. Because is, how- 
ever, relatively small, Patlak's method, for which it is assumed A;4 = 0, is 
still a suitable graphical method for quantifying the uptake rate. 
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Table 1: Brain regions and rate constantsT^ for slice 64 of the Zubal head phantom [j^ . 



Brain Regions 


Ki 


k2 


^3 






ml/min/g 


1 / min 


1/min 


1 / min 


frontal lobes 










occipital lobes 










insula cortex 


0.102 


0.130 


0.062 


0.0068 


temporal lobes 










globus pallidus 










thalamus 


0.082 


0.105 


0.060 


0.0068 


putamen 


0.070 


0.070 


0.054 


0.0068 


caudate nucleus 










internal capsule 










corpus coUosum 


0.054 


0.109 


0.045 


0.0058 


other white matter 











The noise-free input function, with values given in kBq/ml, is given by 
the formulation introduced in 

tG [0,0.25] 

. 339.03(t - 0.25) t G [0.25,0.4433] 

-214.656t + 160.691 t G [0.4433, 0.65] ^ ' 

21.165e-0-7501(t~0.65)0 2359 ^ ^ Q 

Although any reasonable input function, including clinical plasma samples, 
could be used for the simulation this formulation was validated as providing 
a good approximation to the plasma samples of a healthy subject. Given 
the input function, the exact phantom can then be generated using the rate 



constants for the structures detailed in table [H [23j]. The output time frames 
were generated assuming time frames with durations Atj, j = I,-- - ,m, 
m = 22, given in minutes, 0.2, 8 x 0.0333, 2 x 0.1667, 0.2, 0.5, 2 x 1, 2 x 1.5, 
3.5, 2 X 5, 10 and 30. 

In our experiments, in order to control the computational overhead of the 
reconstructions, we reduced the image size of the phantom from 256 x 256 to 
128 X 128. To do this we needed to relabel the structure assigned to a given 
voxel. This was done by averaging the quantity Kik^/{k2 + k^) over the 2- 
by-2 neighboring voxels at the finer resolution. Then the voxel at the coarser 
level was labeled as belonging to the structure for which this average is clos- 
est to the structure value. The kinetic values were then assigned, the TTAC 
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output values calculated, and then projected to the sinogram space yielding 
projected noise-free sinogram data P. For simplification, the instrumental 
and physical effects, including attenuation, Compton scattering, decay and 
random coincidences, were not simulated. Poisson noise was then added to 



the projected data using S = poissrnd(P) where the Matlab@,[19(] func- 
tion poissrnd uses vector P as the means of Poisson densities to generate 
the noisy sinogram data S. Based on this noisy sinogram data, concen- 
tration images are reconstructed using the Expectation-Maximization (EM) 



algorithm, [14 1 



Several data sets were generated: To investigate errors introduced by 
violation of the irreversibility assumption, ^4 = 0, we tested both /C4 = and 
ki > 0. Similarly, we tested both with and without noise in the sinograms to 
investigate the effects of the proposed method due to noise in the sinograms. 
Gaussian noise with noise levels 0%, 5%, 10%, 15% and 20% was added to 
the input function, i.e. 

Cp{t,) = u{tj){l + CYvj), (11) 

where rjj is selected from a standard normal distribution (G(0, 1)), and CV = 
0, 0.05, 0.10, 0.15 and 0.20. 100 random realizations are tested in each case. 

We summarize the test data sets in table O using a character triple 
to classify each test. The first character of the triple indicates whether 
A;4 = or /c4 > 0, or +, respectively. The second character indicates 
the noise level on the input function, or -|-, for noise-free, or with added 
noise, respectively. To indicate the noise level -|- is replaced by 1, 2, 3 and 
4 to indicate noise levels 5%, 10%, 15% and 20% as necessary. The third 
character indicates whether noise is added to the sinogram, again or +, 
respectively. Therefore, for example, +3+ represents the test case for > 0, 
15% noise in u{t) and Poisson noise added to the sinogram. 

4. Results 

To evaluate the simulations quantitatively we define the relative error of 
the k^^ realization for voxel i: 

nk= ' i = l,--- ,iV, A: = l,--- ,100, 

true 

where K^'^ is the estimated value of the true value of the uptake rate i^true 
the k^^ realization for voxel i. Over 100 random simulations, and voxels, 
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Table 2: Summary of the test cases. Here in columns two and three indicates the 
noise-free case, while + indicates noise was added. In column one the indicates ki — 0. 
The same comments apply, but with the irreversibility assumption violated for all triples 
starting with +, indicating k4 > 0. 



^4 


u{t) 


Sinogram 


Comments 










Errors only caused by reconstruction 








+ 


Errors in TTACs caused by reconstruction plus Poisson noise 




+ 





Errors in Cp{t) 




+ 


+ 


General case for irreversible compartmental model 



we calculate the bias, i.e. the mean of the relative errors, and the deviation 
from the mean: 

imN ' looiv ■ ^ ' 

Absolute relative errors \rik\ and associated mean R = Yld=i Tl/k=i kifcl/lOO-^i 
and deviation, D = Yl^^i Ylrk=i Wik\ — ^ are also calculated. Note the 
deviations are li measurements which do not overweight outlier and large 
error samples, as is the case for the /2-based measurements such as the root 
mean squared error. 

In the images shown in the figures we illustrate the calculated uptake 
rates K of the FDG. Images for the CMRGlc can be obtained by directly 
scaling K. In figure [T] we compare the result of using Patlak and TV-Patlak 
for estimating the uptake rates with respect to no noise, 20% noise in the 
input function, Poisson noise in the sinogram, and finally with respect to the 
case in which the irreversibility assumption is violated but without noise in 
the sinogram or input data. In each case the histogram of the relative errors 
is given on the left, the Patlak image in the middle and the TV-Patlak 
on the right. The different scales in the histograms are due to the total 
number of results illustrated. When there is no noise (triples 000 and +00) 
the histogram illustrates results over all voxels but only one simulation, 
while for the noisy simulations the results are for all voxels over all 100 
realizations of the noise. The TV-Patlak images are more homogeneous in 
all cases and the relative errors are smaller. The figures clearly show the 
improvements of employing the TV-Patlak method as compared to using 
Patlak independently for each voxel. This is confirmed in figure [2] in which 
images with noise in the sinogram, positive and different noise levels in 
the input function are shown. 
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Figure 1: Comparison of Patlak and TV-Patlak for imaging the uptake rate K. In each 
case on the left is the histogram for the relative error as compared to the exact simulated 
value, in the middle the Patlak image and on the right the TV-Patlak image. The first 
row compares the estimation by Patlak and TV-Patlak with no noise added. The second, 
third and fourth rows provide the comparison with 20% noise added in the input, Poisson 
noise added to the sinogram, and positive ki, respectively. 



Figure 2: Comparison of Patlak and TV-Patlak for imaging the uptake rate K. In each 
case on the left is the histogram for the relative error as compared to the exact simulated 
value, in the middle the Patlak image and on the right the TV-Patlak image. In each case 
Poisson noise is added to the sinogram and ki > 0. The level of noise for 100 realizations 
added to the input is 5%, 10% and 20% respectively, for each row from top to bottom. 
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Figure 3: Comparison of Patlak (upper-left), TV-Patlak (upper-right), Patlak-GF 
(bottom-left) and Patlak-MF (bottom-right) for imaging the uptake rate K for one sim- 
ulated data case +4-I-, i.e. ki > 0, 20% noise in the input function and Poisson noise in 
the sinogram. 




Relative error r 

|| 

Figure 4; Histograms for the relative error (|12|l in the uptake rate calculated by Patlak, 
TV-Patlak, Patlak-GF and Patlak-MF, for 100 realizations of the data case 4-4-1-, i.e. 
ki > 0, 20% noise in the input function and Poisson noise in the sinogram. 
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Quantitative measurements, confirming the illustrations, are presented 
in table [3l There we also present the results for conventional Patlak's method 
with post-smoothing by two standard filters: 

1. a Gaussian filter (Patlak-GF), size 3-by-3 with standard deviation 0.5, 
generated by Matlab functions f ilter2(f special('gaussian', 3, 0.5), 
img) and 

2. a median filter (Patlak-MF) generated by Matlab function medfi It 2 (img) 
for a 3-by-3 neighborhood. 



Consistent with the observation in 2l|, ll6|] , we find that violation of the Pat- 



lak assumption, = 0, introduces about 10% bias; r f« when ^4 = but 
f « 10% for > 0. The rows TV (std) and "# 10% (#15%)" provide com- 
plementary supporting information, indicating that the TV is minimized by 
TV-Patlak; as compared to Patlak, Patlak-GF and Patlak-MF the number 
of voxels with larger error is reduced. In particular, we emphasize that TV- 
Patlak provides a better noise removal mechanism than popular 
post-filtering approaches. 

In figures [3] and H] we illustrate the uptake rates and relative error in the 
uptake rates, respectively, calculated by Patlak, TV-Patlak, Patlak-GF and 
Patlak-MF for one simulated data case +4+, i.e. k^ > 0, 20% noise in the 
input function and Poisson noise in the sinograms. The uptake rate image 
generated by Patlak-MF is visually smoother than that by TV-Patlak, but 
the equivalent histograms show that the relative error is higher for Patlak- 
MF than for TV-Patlak; the Patlak-MF image is over-smoothed. 

In figure [5] we illustrate the relative error of noisy case +4+ for gray 
matter regions in the phantom, including frontal lobes, occipital lobes, in- 
sula cortex, temporal lobes and globus pallidus, which are of interest for 
Alzheimer's disease research. The error bars show that all estimation meth- 
ods for gray matter regions have negative bias and there are fewer cases with 
high error by TV-Patlak. 

Finally, we note that the computational cost for the TV-Patlak Algo- 
rithm [1] is about 6.5 seconds while for Patlak, Patlak-GF and Patlak-MF 
they are 0.48, 0.55 and 0.58 seconds respectively for a 2D image on a PC 
with IGHz CPU and Matlab code. 

5. Discussion 

In this section we discuss issues relevant to the proposed TV-Patlak 
method. 
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Table 3: Results of phantom simulation: comparison of relative errors between the Patlak 
(P), TV-Patlak (TV-P), Patlak-GF (P-GF) and Patlak- MF (P-MF) methods. TV is the 
mean of the total variation over 100 simulations in each case, std denotes its standard 
deviation. The rows " #10%(#15%)" report the average number of voxels with high 
errors, \rik \ > 10% (\rik \ > 15%) per image. The first row indicates the specific simulation. 



Data case 


000 


040 


00+ 


+00 


+1+ 


+2+ 


+4+ 




r 


-0.009 


0.005 


-0.009 


-0.105 


-0.105 


-0.103 


-0.088 




d 


0.116 


0.136 


0.157 


0.106 


0.148 


0.151 


0.165 




B. 


0.116 


0.136 


0.157 


0.147 


0.176 


0.177 


0.182 


P 


D 


0.088 


0.094 


0.101 


0.081 


0.104 


0.106 


0.113 




TV 


19.488 


22.816 


24.811 


17.246 


22.901 


23.526 


26.074 




std 


0.000 


0.193 


0.235 


0.000 


0.219 


0.225 


0.301 




#10% 


1429 


1718 


2020 


2099 


2272 


2265 


2238 




#15% 


995 


1166 


1453 


1285 


1733 


1732 


1718 




f 


-0.005 


0.020 


0.001 


-0.101 


-0.094 


-0.088 


-0.063 




d 


0.104 


0.120 


0.128 


0.097 


0.121 


0.123 


0.133 


T 


B. 


0.105 


0.121 


0.128 


0.145 


0.157 


0.155 


0.148 


V 


D 


0.095 


0.096 


0.099 


0.078 


0.089 


0.091 


0.096 


1 


TV 


15.663 


17.624 


15.484 


13.650 


13.607 


13.988 


15.344 


P 


std 


0.000 


0.156 


0.220 


0.000 


0.241 


0.254 


0.251 




#10% 


1345 


1462 


1522 


1905 


2164 


2097 


1923 




#15% 


974 


1046 


1102 


1163 


1441 


1415 


1348 




f 


-0.012 


0.002 


-0.012 


-0.108 


-0.108 


-0.106 


-0.091 




d 


0.124 


0.137 


0.151 


0.115 


0.142 


0.143 


0.151 


P 


B 


0.125 


0.137 


0.152 


0.162 


0.175 


0.175 


0.173 


1 


D 


0.106 


0.106 


0.104 


0.088 


0.101 


0.102 


0.106 


G 


TV 


19.052 


20.635 


22.254 


16.881 


20.231 


20.502 


21.653 


F 


std 


0.000 


0.126 


0.171 


0.000 


0.147 


0.155 


0.201 




#10% 


1492 


1589 


1874 


2240 


2297 


2279 


2205 




#15% 


1183 


1191 


1342 


1328 


1722 


1713 


1658 




f 


-0.017 


-0.002 


-0.021 


-0.112 


-0.116 


-0.113 


-0.095 




d 


0.125 


0.138 


0.148 


0.117 


0.139 


0.140 


0.147 


P 


B^ 


0.128 


0.138 


0.150 


0.169 


0.182 


0.181 


0.176 


1 


D 


0.113 


0.115 


0.114 


0.096 


0.103 


0.104 


0.109 


M 


TV 


15.938 


16.295 


16.795 


14.104 


15.026 


15.121 


15.506 


F 


std 


0.000 


0.119 


0.130 


0.000 


0.137 


0.143 


0.166 




#10% 


1499 


1522 


1709 


2342 


2396 


2362 


2215 




#15% 


1118 


1152 


1251 


1335 


1730 


1710 


1615 
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Regularization constant a: There are many approaches for determin- 
ing the choice of an ap propriate regularization constant a. A good 
reference would be [24i | in which the methods of unbiased predictive 
risk-estimation, generalized cross validation, and the L-curve are de- 
scribed. In general, the choice of a depends on the noise level of the 
problem. Here a balances the homogeneity of the uptake rate against 
the residual of the traditional Patlak least squares data fit. For the 
PET imaging application, noise in the TTAC data depends on the 
scanner, the reconstruction method, the tracer dosage and even the 
kinetics of individuals. It is therefore possible to make a standard 
parameter setting for commonly-used environments. The most conve- 



nient method for the selection of a is the so-called L-curve, [10|, 111! ], 
which plots total variation against J^^i \\Wi{A{xi,yi)'^ — bj)||| for all 
tested a. The L-curve clearly displays the compromise between of the 
homogeneity and the residual of the Patlak fitting equations. The a 
corresponding to the left lower corner of the L-curve is considered as 
the optimal choice. One representative L-curve of our simulations is 
illustrated in figure [H We found for our simulations that a suitable 
choice is a = 0.2, but certainly it will in general depend on the recon- 
struction algorithm. For example, the simple EM method and filtered 
backprojection algorithms introduce more noise than the maximum a 
posteriori (MAP) algorithm, [l[ . For real data not only are there 
additional sources of noise but the choice for a will also depend on the 
specific tracer. However, once an appropriate a is found by L-curve 
for a specific imaging environment, it can be fixed for future imaging 
calculations. 

2. Constant f3: We use Newton's method for the case /? 7^ to solve the 
optimization problem ([5]). Other methods for solving the TV prob- 
lem with /3 7^ are discussed in [2j], which includes the primal-dual 
method [J|. A good choice of /3 avoids numerical difficulties for small 
derivatives and provides a good approximation of the TV. For our 
tests we found that the results are not sensitive to the choice of (3 and 
thus suggest the use of /3 = 10~^ as a good choice for FDG-PET brain 
imaging. If we wish to avoid the choice for f3 the TV term in ([5]) can 
be reformulated as a set of linear constraints, and other algorithms are 
possible, 0]. 

3. Boundary values: At each iteration of the algorithm, after updating 
X and before doing any function evaluation or other calculation, the 
boundary voxels need to be updated. The value of each boundary voxel 
is set to the average of its active neighbors in four directions. 
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4. Computational efficiency: For computational efficiency and to min- 
imize memory usage, it is important to not only use sparse storage 
strategies for the relevant matrices but also to use appropriate algo- 
rithms for the solution of the large scale sparse linear systems given 
by i^. Here we use the QMR iteration, 0]. Moreover, if a gen- 
eral Newton's or quasi-Newton method with approximated Hessian 
were used, the cost in time and memory would be much more ex- 
pensive because the approximated Hessian matrix is generally dense. 
In addition to achieving sparsity, the Hessian matrix is calculated 
accurately; and the convergence should be faster. For our simula- 
tions convergence is achieved in eleven iteration steps on average. 
The Matlab code of the TV-Patlak method can be downloaded from 
http : //math . asu . edu/ "^ongbin 



6. Conclusions 

A qualitative improvement in imaging of PET uptake can be achieved by 
using a global model, with the total variation as a penalty term, to obtain the 
voxel uptake rates. The resultant uptake images have spatial homogeneity 
over brain regions with similar kinetics and distinct edges between brain 
regions that have different kinetics. It is statistically validated that the TV- 
Patlak significantly reduces the relative errors of the calculated parameters 
as compared with those generated by Patlak's graphical method, and post- 
smoothing by Gaussian and median filters. 
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Relative error r. for gay matter regions 

Figure 5: Histograms for the relative error (|12|) in the uptake rate of gray matter regions 
calculated by Patlak, TV-Patlak, Patlak-GF and Patlak-MF, for 100 realizations of the 
data case +4+, i.e. ki > 0, 20% noise in the input function and Poisson noise in the 
sinogram. 
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Figure 6: L-curve for the simulation case '+4+'. Total variation against the sum of 



weighted residuals, i.e. X^^^j^ \\Wi{A{xi,yi) —hi) 
a = 0.2 is associated with the left lower "corner" 



, for a from 0.0498 to 0.8187 is plotted. 
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A. Appendix 

To simplify notation, we assume that the weighting matrix in ^ has 
been absorbed into A and bj.To simpUfy the expressions in the objective 
function ^ we introduce the vector function f (z) = (/i(z), /2(z), • • • , /2Ar(z))^, 
where fi{z) = (/)fc(x)2 for i = 1 . . . iV and fi{z) = \\A[xi_N;yi-N] - ^{-nW^ 
for i = N + 1...2N. 
Gradient calculation 

To find the gradient of ^*(z), we first derive the Jacobian of f(z). 
1. For i = 1, ■ • • , and / = 1, • • • , 2N 



dzi 



4:Xi '^•^if, '^^ir 1 ^ — i 



0, 



if / = ib 
if I = ir 
otherwise. 



(13) 



2. For i = iV + 1, • • • , 2iV and / = 1, • ■ ■ , 2N 

— = { Pi2, it I = t 

0, otherwise, 



9z 



where 



Pil 

Pi2 



2A^A f "^^-^ 

Vi-N 



2A^h. 



■i-N- 



(14) 



Because Vx^i = l/(2(/.,)Vx/i for i = 1, • • • , iV, = (Vf(z))^g, 

where 

g = [l/(2(/>i), 1/(202), ■ ■ ■ , l/(2<^iv), a, a, • • • , af . 
Hessian matrix 

To find the Hessian matrix for *J'(z), we first derive the Hessian matrix for 
each function /j(z). 

1. For i = 1, • • • ,N, Vxx/i is sparse 



fi 



I 


% 






4 


-2 


-2^ 




-2 


2 







-2 





2) 





(15) 
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2. For i = N + 1, ■ ■ ■ , 2N, the Hessian matrix is again sparse. 

i- N i . s 

VL/. = f m m \ ^-N where ( ) = 2^^ A. 

(16) 



Thus 



j=l ^ ^ i=N+l 



where 
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